Fitting models in Pumas
2021-04-20
using Random using Pumas using PlottingUtilities using PumasPlots
Fitting a PK model
In this tutorial we will go through the steps required to fit a model in Pumas.jl.
The dosage regimen
We start by simulating a population from a two compartment model with oral absorption, and then we show how to fit and do some model validation using the fit output.
The dosage regimen is an dose of 100 into Depot at time 0, followed by two additional (addl=2) doses every fourth hour
repeated_dose_regimen = DosageRegimen(100, time = 0, ii = 4, addl = 2)
| time | cmt | amt | evid | ii | addl | rate | duration | ss | route | |
|---|---|---|---|---|---|---|---|---|---|---|
| Float64 | Int64 | Float64 | Int8 | Float64 | Int64 | Float64 | Float64 | Int8 | Route… | |
| 1 | 0.0 | 1 | 100.0 | 1 | 4.0 | 2 | 0.0 | 0.0 | 0 | NullRoute |
As ususal, let's define a function to choose body weight randomly per subject
choose_covariates() = (Wt = rand(55:80),)
choose_covariates (generic function with 1 method)
and generate a population of subjects with a random weight generated from the covariate function above
pop = map(i -> Subject(id = i, events = repeated_dose_regimen, observations = (dv = nothing,), covariates = choose_covariates()), 1:24)
Population Subjects: 24 Covariates: Wt Observables: dv
We now have 24 subjects equipped with a weight and a dosage regimen.
The PK model of drug concentration and elimination
To simulate a data set and attempt to estimate the data generating parameters, we have to set up the actual pharmacokinetics (PK) model and simulate the data. We use the closed form model called Depots1Central1Periph1 which is a two compartment model with first order absorption. This requires CL, Vc, Ka, Vp, and Q to be defined in the @pre-block, since they define the rates of transfer between (and out of) the compartments
mymodel = @model begin @param begin cl ∈ RealDomain(lower = 0.0, init = 1.0) tv ∈ RealDomain(lower = 0.0, init = 10.0) ka ∈ RealDomain(lower = 0.0, init = 1.0) q ∈ RealDomain(lower = 0.0, init = 0.5) Ω ∈ PDiagDomain(init = [0.9,0.07, 0.05]) σ_prop ∈ RealDomain(lower = 0,init = 0.03) end @random begin η ~ MvNormal(Ω) end @covariates Wt @pre begin CL = cl * (Wt/70)^0.75 * exp(η[1]) Vc = tv * (Wt/70) * exp(η[2]) Ka = ka * exp(η[3]) Vp = 30.0 Q = q end @dynamics Depots1Central1Periph1 @derived begin cp := @. 1000*(Central / Vc) # We use := because we don't want simobs to store the variable dv ~ @. Normal(cp, abs(cp)*σ_prop) end end
PumasModel Parameters: cl, tv, ka, q, Ω, σ_prop Random effects: η Covariates: Wt Dynamical variables: Depot, Central, Peripheral Derived: dv Observed: dv
Some parameters are left free by giving them domains in the @param-block, and one PK parameter (the volume of distribution of the peripheral compartment) is fixed to 20.0.
Simulating the individual observations
The simobs function is used to simulate individual time series. We input the model, the population of Subjects that currently only have dosage regimens and covariates, the parameter vector and the times where we want to simulate. Since we have a proportional error model we avoid observations at time zero to avoid degenerate distributions of the dependent variable. The problem is, that if the concentration is zero the variance in distribution of the explained variable will also be zero. Let's use the default parameters, as set in the @param-block, and simulate the data
param = init_params(mymodel) Random.seed!(1234) obs = simobs(mymodel, pop, param, obstimes = 1:1:72)
interactive!(false) p = sim_plot(mymodel, obs, observations = :dv) p[1]
Fitting the model
To fit the model, we use the fit function. It requires a model, a population, a named tuple of parameters and a likelihood approximation method.
result = fit(mymodel, Subject.(obs), param, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter Function value Gradient norm
0 1.002801e+04 1.608331e+02
* time: 6.198883056640625e-5
1 1.002794e+04 4.438569e+01
* time: 0.055644989013671875
2 1.002778e+04 3.281777e+01
* time: 0.06789994239807129
3 1.002723e+04 1.152964e+01
* time: 0.07953786849975586
4 1.002646e+04 1.498682e+01
* time: 0.0889890193939209
5 1.002549e+04 3.528519e+01
* time: 0.0987548828125
6 1.002525e+04 5.396675e+00
* time: 0.10743594169616699
7 1.002515e+04 1.900268e+00
* time: 0.1170189380645752
8 1.002507e+04 4.798916e+00
* time: 0.12464499473571777
9 1.002495e+04 8.321538e+00
* time: 0.13320684432983398
10 1.002487e+04 5.994750e+00
* time: 0.14234185218811035
11 1.002485e+04 1.719327e+00
* time: 0.15235590934753418
12 1.002485e+04 1.420724e-01
* time: 0.15978598594665527
13 1.002485e+04 6.651317e-02
* time: 0.16875195503234863
14 1.002485e+04 2.056923e-01
* time: 0.17668485641479492
15 1.002485e+04 2.617325e-01
* time: 0.18510890007019043
16 1.002485e+04 1.458920e-01
* time: 0.19348788261413574
17 1.002485e+04 3.011164e-02
* time: 0.20208191871643066
18 1.002485e+04 1.971643e-03
* time: 0.2097618579864502
19 1.002485e+04 2.117735e-03
* time: 0.21756601333618164
20 1.002485e+04 4.095262e-03
* time: 0.2253720760345459
21 1.002485e+04 4.289614e-03
* time: 0.23316693305969238
22 1.002485e+04 2.136597e-03
* time: 0.24097490310668945
23 1.002485e+04 1.319400e-03
* time: 0.2505300045013428
24 1.002485e+04 1.183455e-03
* time: 0.2590649127960205
25 1.002485e+04 1.183389e-03
* time: 0.2725410461425781
26 1.002485e+04 9.994834e-04
* time: 0.28187084197998047
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10024.853
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.7039
tv 9.774
ka 1.0298
q 0.50035
Ω₁,₁ 0.76081
Ω₂,₂ 0.048151
Ω₃,₃ 0.058265
σ_prop 0.03027
---------------------
Of course, we started the fitting at the true parameters, so let us define our own starting parameters, and fit based on those values
alternative_param = ( cl = 0.5, tv = 9.0, ka = 1.3, q = 0.3, Ω = Diagonal([0.18,0.04, 0.03]), σ_prop = 0.04) fit(mymodel, Subject.(obs), alternative_param, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter Function value Gradient norm
0 2.192254e+04 3.944925e+04
* time: 4.982948303222656e-5
1 1.396191e+04 1.736510e+04
* time: 0.03451800346374512
2 1.106810e+04 4.365316e+03
* time: 0.04655599594116211
3 1.078869e+04 1.310417e+03
* time: 0.05642580986022949
4 1.072387e+04 1.278217e+03
* time: 0.06484198570251465
5 1.065219e+04 1.840211e+03
* time: 0.07382488250732422
6 1.041598e+04 5.547035e+03
* time: 0.08307480812072754
7 1.009012e+04 3.261406e+03
* time: 0.09153485298156738
8 1.006355e+04 3.884887e+02
* time: 0.10253190994262695
9 1.006249e+04 1.069016e+02
* time: 0.11028385162353516
10 1.006200e+04 1.048583e+02
* time: 0.11846089363098145
11 1.006020e+04 2.220255e+02
* time: 0.12783384323120117
12 1.005671e+04 3.359330e+02
* time: 0.136185884475708
13 1.005210e+04 3.203693e+02
* time: 0.1452028751373291
14 1.004922e+04 1.521796e+02
* time: 0.1540839672088623
15 1.004858e+04 3.634242e+01
* time: 0.16219496726989746
16 1.004844e+04 3.631365e+01
* time: 0.1716480255126953
17 1.004830e+04 4.239529e+01
* time: 0.17960286140441895
18 1.004791e+04 9.158692e+01
* time: 0.187391996383667
19 1.004701e+04 1.583750e+02
* time: 0.19507694244384766
20 1.004501e+04 2.382150e+02
* time: 0.2052619457244873
21 1.004146e+04 2.906203e+02
* time: 0.21512103080749512
22 1.003672e+04 2.492799e+02
* time: 0.22508692741394043
23 1.003337e+04 1.015827e+02
* time: 0.23441290855407715
24 1.003270e+04 2.552435e+01
* time: 0.24299883842468262
25 1.003265e+04 2.654450e+01
* time: 0.25171995162963867
26 1.003259e+04 2.677991e+01
* time: 0.26029491424560547
27 1.003242e+04 2.573818e+01
* time: 0.2692899703979492
28 1.003209e+04 3.439711e+01
* time: 0.27742981910705566
29 1.003152e+04 4.337631e+01
* time: 0.28592586517333984
30 1.003090e+04 3.629501e+01
* time: 0.29424381256103516
31 1.003058e+04 1.722268e+01
* time: 0.30318784713745117
32 1.003049e+04 1.913596e+01
* time: 0.31111884117126465
33 1.003046e+04 1.829670e+01
* time: 0.31977081298828125
34 1.003041e+04 1.664881e+01
* time: 0.3277459144592285
35 1.003029e+04 1.962906e+01
* time: 0.33602285385131836
36 1.003001e+04 3.529991e+01
* time: 0.3442549705505371
37 1.002941e+04 5.529316e+01
* time: 0.3540019989013672
38 1.002836e+04 7.062268e+01
* time: 0.3624589443206787
39 1.002697e+04 6.460823e+01
* time: 0.3717219829559326
40 1.002578e+04 3.541353e+01
* time: 0.3831770420074463
41 1.002521e+04 9.253277e+00
* time: 0.39461803436279297
42 1.002508e+04 3.191871e+00
* time: 0.406674861907959
43 1.002505e+04 2.254884e+00
* time: 0.41837596893310547
44 1.002505e+04 1.944729e+00
* time: 0.431121826171875
45 1.002505e+04 1.945441e+00
* time: 0.44263291358947754
46 1.002505e+04 1.939199e+00
* time: 0.45410799980163574
47 1.002504e+04 1.925990e+00
* time: 0.46538400650024414
48 1.002504e+04 1.906824e+00
* time: 0.4767138957977295
49 1.002504e+04 1.868390e+00
* time: 0.4896550178527832
50 1.002504e+04 1.794271e+00
* time: 0.501863956451416
51 1.002503e+04 1.737884e+00
* time: 0.5245168209075928
52 1.002501e+04 2.818073e+00
* time: 0.5375969409942627
53 1.002496e+04 3.787146e+00
* time: 0.5500078201293945
54 1.002490e+04 3.550240e+00
* time: 0.562675952911377
55 1.002486e+04 1.667143e+00
* time: 0.5786800384521484
56 1.002485e+04 3.962635e-01
* time: 0.5912308692932129
57 1.002485e+04 7.931606e-02
* time: 0.6037249565124512
58 1.002485e+04 3.453083e-02
* time: 0.6149590015411377
59 1.002485e+04 2.767901e-02
* time: 0.6278948783874512
60 1.002485e+04 2.767901e-02
* time: 0.6490998268127441
61 1.002485e+04 2.767901e-02
* time: 0.6750028133392334
62 1.002485e+04 2.769532e-02
* time: 0.6898720264434814
63 1.002485e+04 2.769532e-02
* time: 0.717491865158081
64 1.002485e+04 1.259742e-01
* time: 0.7309579849243164
65 1.002485e+04 6.421964e-02
* time: 0.7441809177398682
66 1.002485e+04 3.399018e-02
* time: 0.7608239650726318
67 1.002485e+04 1.673434e-01
* time: 0.7866170406341553
68 1.002485e+04 3.342486e-01
* time: 0.8018498420715332
69 1.002485e+04 4.374753e-01
* time: 0.8140339851379395
70 1.002485e+04 3.371089e-01
* time: 1.5997569561004639
71 1.002485e+04 1.180698e-01
* time: 1.6088249683380127
72 1.002485e+04 1.079942e-02
* time: 1.6190450191497803
73 1.002485e+04 1.079942e-02
* time: 1.6366229057312012
74 1.002485e+04 1.079942e-02
* time: 1.656303882598877
75 1.002485e+04 1.079906e-02
* time: 1.6758918762207031
76 1.002485e+04 1.079905e-02
* time: 1.6911888122558594
77 1.002485e+04 1.079901e-02
* time: 1.7079260349273682
78 1.002485e+04 1.079897e-02
* time: 1.7238948345184326
79 1.002485e+04 1.079896e-02
* time: 1.7403960227966309
80 1.002485e+04 1.079896e-02
* time: 1.7574079036712646
81 1.002485e+04 1.079895e-02
* time: 1.7761459350585938
82 1.002485e+04 1.079895e-02
* time: 1.7963709831237793
83 1.002485e+04 1.079895e-02
* time: 1.8134968280792236
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10024.853
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 0.7039
tv 9.774
ka 1.0298
q 0.50035
Ω₁,₁ 0.76079
Ω₂,₂ 0.048149
Ω₃,₃ 0.058264
σ_prop 0.03027
---------------------
and we see that the estimates are essentially the same up to numerical noise.
To augment the basic information listed when we print the results, we can use infer to provide RSEs and confidence intervals
infer(result)
Asymptotic inference results using sandwich estimator
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -10024.853
Number of subjects: 24
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 0.7039 0.12541 [0.45809 ; 0.9497 ]
tv 9.774 0.43837 [8.9148 ; 10.633 ]
ka 1.0298 0.052264 [0.92735 ; 1.1322 ]
q 0.50035 0.0011047 [0.49818 ; 0.50251 ]
Ω₁,₁ 0.76081 0.2891 [0.19418 ; 1.3274 ]
Ω₂,₂ 0.048151 0.013182 [0.022315; 0.073986]
Ω₃,₃ 0.058265 0.016079 [0.02675 ; 0.08978 ]
σ_prop 0.03027 0.00051989 [0.029251; 0.031289]
--------------------------------------------------------------------
So as we observed earlier, the parameters look like they have sensible values. The confidence intervals are a bit wide, and especially so for the random effect variability parameters. To see how we can use simulation to better understand the statistical properties of our model, we can simulate a much larger population and try again
pop_big = map(i -> Subject(id = i, events = repeated_dose_regimen, observations = (dv = nothing,), covariates = choose_covariates()), 1:100) obs_big = simobs(mymodel, pop_big, param, obstimes=1:1:72) result_big = fit(mymodel, Subject.(obs_big), param, Pumas.FOCEI(), ensemblealg = EnsembleThreads()) infer(result_big)
Iter Function value Gradient norm
0 3.679364e+04 7.630791e+02
* time: 8.20159912109375e-5
1 3.679332e+04 1.425848e+02
* time: 0.16926002502441406
2 3.679307e+04 1.234952e+02
* time: 0.21648120880126953
3 3.679285e+04 5.195340e+01
* time: 0.7413251399993896
4 3.679173e+04 7.220420e+01
* time: 0.7796931266784668
5 3.679106e+04 4.575628e+01
* time: 0.8308331966400146
6 3.679072e+04 1.144429e+01
* time: 0.869042158126831
7 3.679068e+04 4.149921e+00
* time: 0.9115841388702393
8 3.679067e+04 5.467597e-01
* time: 0.9500491619110107
9 3.679067e+04 5.097025e-01
* time: 0.9821732044219971
10 3.679066e+04 7.409105e-01
* time: 1.0126211643218994
11 3.679066e+04 5.546808e-01
* time: 1.0436930656433105
12 3.679066e+04 1.370981e-01
* time: 1.0739221572875977
13 3.679066e+04 1.578532e-02
* time: 1.1063601970672607
14 3.679066e+04 1.578538e-02
* time: 1.1588921546936035
15 3.679066e+04 1.581460e-02
* time: 1.2266860008239746
16 3.679066e+04 1.203715e-01
* time: 1.2714600563049316
17 3.679066e+04 1.505423e-01
* time: 1.3377540111541748
Asymptotic inference results using sandwich estimator
Successful minimization: false
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -36790.662
Number of subjects: 100
Number of parameters: Fixed Optimized
0 8
Observation records: Active Missing
dv: 7200 0
Total: 7200 0
--------------------------------------------------------------------
Estimate SE 95.0% C.I.
--------------------------------------------------------------------
cl 1.215 0.11821 [0.98329 ; 1.4467 ]
tv 10.041 0.27033 [9.5115 ; 10.571 ]
ka 1.0159 0.023039 [0.9707 ; 1.061 ]
q 0.5004 0.00051675 [0.49939 ; 0.50141 ]
Ω₁,₁ 0.94724 0.13267 [0.68721 ; 1.2073 ]
Ω₂,₂ 0.071929 0.010663 [0.051029; 0.092829]
Ω₃,₃ 0.049596 0.006597 [0.036666; 0.062526]
σ_prop 0.029833 0.00024886 [0.029345; 0.030321]
--------------------------------------------------------------------
This time we see similar estimates, but much narrower confidence intervals across the board.
Estimating a misspecified model
To explore some of the diagnostics tools available in Pumas, we can try to set up a model that does not fit out data generating process. This time we propose a one compartent model. The problem with estimating a one compartment model when the data comes from a two compartment model, is that we cannot capture the change in slope on the concentration profile you get with a two compartment model. This means that even if we can capture the model fit someone well on average, we should expect to see systematic trends in the residual diagnostics post estimation.
mymodel_misspec = @model begin @param begin cl ∈ RealDomain(lower = 0.0, init = 1.0) tv ∈ RealDomain(lower = 0.0, init = 20.0) ka ∈ RealDomain(lower = 0.0, init = 1.0) Ω ∈ PDiagDomain(init = [0.12, 0.05, 0.08]) σ_prop ∈ RealDomain(lower = 0, init = 0.03) end @random begin η ~ MvNormal(Ω) end @pre begin CL = cl * (Wt/70)^0.75 * exp(η[1]) Vc = tv * (Wt/70) * exp(η[2]) Ka = ka * exp(η[3]) end @covariates Wt @dynamics Depots1Central1 @derived begin cp = @. 1000*(Central / Vc) dv ~ @. Normal(cp, abs(cp)*σ_prop) end end
PumasModel Parameters: cl, tv, ka, Ω, σ_prop Random effects: η Covariates: Wt Dynamical variables: Depot, Central Derived: cp, dv Observed: cp, dv
alternative_param_no_q = ( cl = 0.5, tv = 9.0, ka = 1.3, Ω = Diagonal([0.18,0.04, 0.03]), σ_prop = 0.04) result_misspec = fit(mymodel_misspec, Subject.(obs), alternative_param_no_q, Pumas.FOCEI(), ensemblealg = EnsembleThreads())
Iter Function value Gradient norm
0 9.617797e+04 1.693060e+05
* time: 8.0108642578125e-5
1 2.377470e+04 2.167379e+04
* time: 0.040180206298828125
2 2.104349e+04 1.578102e+04
* time: 0.0628671646118164
3 1.691859e+04 6.360556e+03
* time: 0.0836172103881836
4 1.568972e+04 3.107010e+03
* time: 0.10488510131835938
5 1.514405e+04 1.262513e+03
* time: 0.12670207023620605
6 1.499042e+04 4.451224e+02
* time: 0.1480541229248047
7 1.495861e+04 4.080309e+02
* time: 0.16955018043518066
8 1.495423e+04 4.069141e+02
* time: 0.19259309768676758
9 1.495225e+04 4.046033e+02
* time: 0.2149050235748291
10 1.494680e+04 3.956093e+02
* time: 0.23790311813354492
11 1.493489e+04 3.723547e+02
* time: 0.25963521003723145
12 1.490880e+04 3.650006e+02
* time: 0.2830040454864502
13 1.486507e+04 4.361740e+02
* time: 0.3068711757659912
14 1.481223e+04 4.027373e+02
* time: 0.330794095993042
15 1.477083e+04 2.426534e+02
* time: 0.3521871566772461
16 1.475778e+04 7.890989e+01
* time: 0.37361812591552734
17 1.475683e+04 8.094033e+01
* time: 0.3956911563873291
18 1.475655e+04 8.196633e+01
* time: 0.4183080196380615
19 1.475644e+04 8.206960e+01
* time: 0.44081997871398926
20 1.475593e+04 8.192692e+01
* time: 0.4632279872894287
21 1.475485e+04 8.077322e+01
* time: 0.4852120876312256
22 1.475195e+04 7.638707e+01
* time: 0.509768009185791
23 1.474556e+04 6.545525e+01
* time: 0.5339281558990479
24 1.473330e+04 7.169812e+01
* time: 0.5564100742340088
25 1.471788e+04 6.294130e+01
* time: 0.5815870761871338
26 1.470981e+04 2.802554e+01
* time: 0.6071221828460693
27 1.470869e+04 1.324312e+01
* time: 0.6262261867523193
28 1.470863e+04 1.357644e+01
* time: 0.6474640369415283
29 1.470862e+04 1.315681e+01
* time: 0.6711502075195312
30 1.470862e+04 1.270641e+01
* time: 0.6967620849609375
31 1.470862e+04 1.239475e+01
* time: 0.7218751907348633
32 1.470861e+04 1.174411e+01
* time: 0.7490310668945312
33 1.470858e+04 1.098416e+01
* time: 1.4301671981811523
34 1.470851e+04 1.135417e+01
* time: 1.4516420364379883
35 1.470833e+04 1.156782e+01
* time: 1.4727990627288818
36 1.470794e+04 1.097780e+01
* time: 1.4933290481567383
37 1.470725e+04 8.691388e+00
* time: 1.5140891075134277
38 1.470638e+04 6.461460e+00
* time: 1.5351810455322266
39 1.470559e+04 4.954152e+00
* time: 1.556359052658081
40 1.470520e+04 2.765504e+00
* time: 1.5781450271606445
41 1.470508e+04 2.724332e+00
* time: 1.5994701385498047
42 1.470504e+04 2.675423e+00
* time: 1.617218017578125
43 1.470501e+04 2.616206e+00
* time: 1.6358060836791992
44 1.470501e+04 2.563224e+00
* time: 1.6564819812774658
45 1.470500e+04 2.535331e+00
* time: 1.6736981868743896
46 1.470500e+04 2.527151e+00
* time: 1.6900320053100586
47 1.470500e+04 2.521881e+00
* time: 1.7055530548095703
48 1.470500e+04 2.512274e+00
* time: 1.7213730812072754
49 1.470500e+04 2.498372e+00
* time: 1.7394030094146729
50 1.470500e+04 2.473337e+00
* time: 1.756925106048584
51 1.470500e+04 2.426461e+00
* time: 1.7831761837005615
52 1.470498e+04 2.331011e+00
* time: 1.8028221130371094
53 1.470496e+04 2.499398e+00
* time: 1.8230462074279785
54 1.470490e+04 2.693471e+00
* time: 1.8456790447235107
55 1.470480e+04 2.774137e+00
* time: 1.8676741123199463
56 1.470464e+04 2.498396e+00
* time: 1.8870291709899902
57 1.470447e+04 1.656331e+00
* time: 1.908945083618164
58 1.470439e+04 6.916492e-01
* time: 1.9297130107879639
59 1.470438e+04 5.068855e-01
* time: 1.9506402015686035
60 1.470438e+04 5.316519e-01
* time: 1.9687931537628174
61 1.470438e+04 5.312769e-01
* time: 1.9852612018585205
62 1.470438e+04 5.295780e-01
* time: 2.0021719932556152
63 1.470438e+04 5.279847e-01
* time: 2.018195152282715
64 1.470438e+04 5.243690e-01
* time: 2.035902976989746
65 1.470438e+04 5.188250e-01
* time: 2.0537030696868896
66 1.470438e+04 5.088026e-01
* time: 2.0722310543060303
67 1.470438e+04 4.911665e-01
* time: 2.0908701419830322
68 1.470438e+04 7.153115e-01
* time: 2.110969066619873
69 1.470437e+04 1.100789e+00
* time: 2.1326301097869873
70 1.470436e+04 1.597138e+00
* time: 2.1561310291290283
71 1.470434e+04 2.018414e+00
* time: 2.178997039794922
72 1.470431e+04 1.914077e+00
* time: 2.2029471397399902
73 1.470427e+04 9.915496e-01
* time: 2.9200031757354736
74 1.470426e+04 1.919551e-01
* time: 2.9422202110290527
75 1.470425e+04 2.509096e-01
* time: 2.9652931690216064
76 1.470425e+04 3.312967e-01
* time: 2.9842851161956787
77 1.470425e+04 2.571555e-01
* time: 3.003941059112549
78 1.470425e+04 8.701710e-02
* time: 3.025247097015381
79 1.470425e+04 2.446243e-02
* time: 3.0460000038146973
80 1.470425e+04 6.000590e-02
* time: 3.0664451122283936
81 1.470425e+04 5.240235e-02
* time: 3.0837981700897217
82 1.470425e+04 2.019211e-02
* time: 3.101785182952881
83 1.470425e+04 4.286986e-03
* time: 3.1202361583709717
84 1.470425e+04 1.288839e-02
* time: 3.1361382007598877
85 1.470425e+04 1.114816e-02
* time: 3.15051007270813
86 1.470425e+04 4.322042e-03
* time: 3.1664481163024902
87 1.470425e+04 9.483345e-04
* time: 3.1831350326538086
FittedPumasModel
Successful minimization: true
Likelihood approximation: Pumas.FOCE
Log-likelihood value: -14704.247
Number of subjects: 24
Number of parameters: Fixed Optimized
0 7
Observation records: Active Missing
dv: 1728 0
Total: 1728 0
---------------------
Estimate
---------------------
cl 1.0129
tv 21.454
ka 666470.0
Ω₁,₁ 0.22529
Ω₂,₂ 0.03929
Ω₃,₃ 0.053478
σ_prop 0.43937
---------------------
First off, the absorption flow parameters Ka is quite off the charts, so that would be a warning sign off the bat, but let us try to use a tool in the toolbox to asses the goodness of fit. We get these by using the inspect function
res_inspect_2cmp = inspect(result) res_inspect_1cmp = inspect(result_misspec) p1 = goodness_of_fit(res_inspect_1cmp) p2 = goodness_of_fit(res_inspect_2cmp)
1-element Vector{AbstractPlotting.Figure}:
Figure()
p1[1]
